#define _USE_MATH_DEFINES
#include <iostream>
#include "matrix.hpp"
#include "matrixdecomposition.hpp"
#include "leastSquareSolver.hpp"
#include "BermudanSwaption.hpp"
#include "swap.hpp"
#include "LMMPricer.hpp"
#include "varianceStructure.hpp"
#include "SimpleMultiDimensionGaussian.hpp"
#include "LMMPathPool.hpp"
#include "LMMBermudanSwaptionPricer.hpp"
#include "BlackScholesPricer.hpp"
#include <fstream>
using namespace std;
using namespace TermStructure;
template<typename T>
std::ostream &outputTo(std::ostream &stream, T &m)
{
	for(int i = 0; i < m.GetLength(0); i++)
	{
		for(int j = 0; j < m.GetLength(1); j++)
		{
			stream<<m(i,j)<<",";
		}
		stream<<std::endl;
	}
	return stream;
}

void t()
{
    
	TimeSequence<double> ts;
	for(int i = 1; i <= 40; i++)
	{
		ts.PutTime(i*0.25);
	}
	ForwardRateSequence<double> frs;
	frs.SetTimeSequence(ts);
	for(int i = 0; i < ts.GetSize(); i++)
	{
		frs[i] = 0.05;
	}
	VarianceStructure<double> vs;
	vs.d=0.05;
	vs.a=0.09;
	vs.b=0.44;
	vs.c=0.11;
	vs.beta2=0.1;
    vs.beta1 = 0;
    Swap<double> swap(ts, 40, 4);
	LMMBermudanSwaptionPricer<VarianceStructure<double>,PCABasedMultiDimensionGaussian> pricer;
	pricer.SetLMMAdvancer(ts, vs);
	for(int i = 0; i < pricer.InitialForwardRateSequence().GetSize(); i++)
	{
		pricer.InitialForwardRateSequence()[i] = 0.05;
	}
	//cout<<pricer.PriceEuropeanSwaption(swap, 0.05, 4);
	LMMPathPool<VarianceStructure<double>,PCABasedMultiDimensionGaussian> lmmpathpool;
	lmmpathpool.SetLMMAdvancer(pricer.UnderlyingAdvancer);
	lmmpathpool.CurrentTimeIndex = 0;
	lmmpathpool.EndIndex = 39;
	lmmpathpool.NumberOfPath = 10000;
	lmmpathpool.Generate();
	lmmpathpool(0,0).EstablishDiscountFactors();
	double t = lmmpathpool(0,0).DiscountFactorAt(4);
	swap.Strike = 0.05;
	double rr = 0;
	for(int i = 0; i < 10000; i++)
	{
		double v = swap.Price(lmmpathpool(i,4));
		if(v >0) rr += v/10000;
	}
	cout<<rr*t;
    
	swap.SwapEnterIndex = 0;
	swap.Strike = 0.05;
	cout<<pricer.PriceBermudanSwaption(swap, 10000);
    BlackScholesPricer<VarianceStructure<double,long>, double> bsp;
    bsp.SetTimeSequence(ts);
    bsp.PriceCaplet(lmmpathpool(0,0), vs, 0.05, 3);
}

int main(int argc, char const * argv[])
{
	TimeSequence<double> ts;
	for(int i = 1; i <= 40; i++)
	{
		ts.PutTime(i*0.25);
	}
	ForwardRateSequence<double> frs;
	frs.SetTimeSequence(ts);
	for(int i = 0; i < ts.GetSize(); i++)
	{
		frs[i] = 0.05;
	}
	VarianceStructure<double> vs;
	vs.d=0.05;
	vs.a=0.09;
	vs.b=0.44;
	vs.c=0.11;
	vs.beta2=0.1;
    vs.beta1 = 0;
	vs.SetTimeSequence(ts);
    BlackScholesPricer<VarianceStructure<double,long>, double> bsp;
    bsp.SetTimeSequence(ts);
    cout<<vs.VolatilityIntegrator(0, 0.75, 4, 4)<<endl;
    cout<<bsp.PriceCaplet(frs, vs, 0.05, 3)<<endl;
	LMMBermudanSwaptionPricer<VarianceStructure<double>,SimpleMultiDimensionGaussian> pricer;
	pricer.SetLMMAdvancer(ts, vs);
	for(int i = 0; i < pricer.InitialForwardRateSequence().GetSize(); i++)
	{
		pricer.InitialForwardRateSequence()[i] = 0.05;
	}
    double r;
    LMMPathPool<VarianceStructure<double>,SimpleMultiDimensionGaussian> lmmpathpool;
	lmmpathpool.SetLMMAdvancer(pricer.UnderlyingAdvancer);
    pricer.UnderlyingAdvancer.AnchorIndex = 40;
	lmmpathpool.CurrentTimeIndex = 0;
	lmmpathpool.EndIndex = 39;
	lmmpathpool.NumberOfPath = 50000;
	lmmpathpool.Generate();
    
    r = 0;
    lmmpathpool(0,0).EstablishDiscountFactors();
    Swap<double> aswap(ts);
    for(int k = 1; k <= 39; k++)
    {
        r = 0;
        aswap.SwapEnterIndex = k; aswap.SwapEndIndex = 40;
        aswap.Strike = 0.06;
    for(int i = 0; i < 50000; i++)
    {
        
        double v = aswap.Price(lmmpathpool(i,k));
        if(v>=0) r += v/50000/lmmpathpool(i,k).DiscountFactorAt(40,k);
    }
        cout<<r*lmmpathpool(0,0).DiscountFactorAt(40)<<endl;
    }
    lmmpathpool(0,0).EstablishDiscountFactors();
    cout<<r*lmmpathpool(0,0).DiscountFactorAt(3)<<endl;
    pricer.UnderlyingAdvancer.SetForwardRateSequence(lmmpathpool(0,0));
    pricer.UnderlyingAdvancer.CurrentTime = 0;
    pricer.UnderlyingAdvancer.AnchorIndex = 3;
    pricer.UnderlyingAdvancer.TimeStep = 0.75;
    pricer.UnderlyingAdvancer.PreAdvance();
    r = 0;
    
    for(int i = 0; i < 100000; i++)
    {
        pricer.UnderlyingAdvancer.SingleAdvance();
        double v = (pricer.UnderlyingAdvancer.NextSequence()[4] - 0.05)*0.25/(1+pricer.UnderlyingAdvancer.NextSequence()[4]*0.25);
        if(v >= 0) r+= v/100000;
    }
    cout<<r*lmmpathpool(0,0).DiscountFactorAt(3)<<endl;
   
    
	lmmpathpool(0,0).EstablishDiscountFactors();
	double t = lmmpathpool(0,0).DiscountFactorAt(4);
    //cout<<pricer.PriceCaplet(0, 0.75, 1, 0.05,1000000);
	return 0;
}